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Abstract: We advocate a method to improve systematically the self-consistent harmonic 
approximation (or the Gaussian approximation), which has been employed extensively in 
condensed matter physics and statistical mechanics. Such a method was previously applied 
to the IIB matrix model, a conjectured nonperturbative definition of type IIB superstring 
theory in ten dimensions. Remarkably the dominance of four-dimensional space-time in 
the partition function was suggested from calculations up to the 3rd order. Recently this 
calculation has been extended to the 5th order, and the same conclusion has been obtained. 
Here we apply this Gaussian expansion method to the bosonic version of the IIB matrix 
model, where Monte Carlo results are available, and demonstrate the convergence of the 
method by explicit calculations up to the 7th order. More generally we study matrix models 
obtained from dimensional reduction of SU(iV) Yang-Mills theory in D dimensions, where 
the D = 10 case corresponds to the bosonic IIB matrix model. Convergence becomes faster 
as D increases, and for D > 10 it is already achieved at the 3rd order. 
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1. Introduction 

It is well appreciated that perturbation theory has a limited range of applicability, although 
its importance in theoretical physics can never be overstated. First of all, perturbative 
expansions in most cases yield merely an asymptotic series, which starts to diverge at some 
finite order. There are also situations in which the expansion parameter is too large to 
make perturbative calculations reliable or even meaningful. 

Superstring theory provides such a situation where nonperturbative effects are con- 
sidered to be extremely important. From this point of view, the recent proposals for its 
nonperturbative formulations should be considered as a substantial progress. For instance 
the IIB matrix model [1] is conjectured to be a nonperturbative definition of type IIB 
superstring theory in ten dimensions, and it is speculated to explain the dynamical origin 
of the Standard Model including the space-time dimensionality, the gauge group and the 
number of generations. Much effort has been made to understand the dynamics of this 
model. In particular, Monte Carlo simulations may be as useful as in lattice gauge theo- 
ries, and indeed various simplified versions of the IIB matrix model have been successfully 
studied [2, 3]. However, it is not straightforward to extend these studies to the IIB matrix 
model, where the complex fermion determinant causes the notorious complex-action prob- 
lem. In fact it is speculated that the phase of the fermion determinant plays a crucial role in 
the dynamical generation of four-dimensional space-time [4]. Recently a new Monte Carlo 
technique [5] is developed to include the effect of the phase. Preliminary results concerning 
the space-time dimensionality are encouraging, but it remains to be seen whether definite 
conclusions can be reached along this line. In any case it would be certainly desirable to 
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develop an alternative method which enables a nonperturbative access to the dynamics of 
the model. 

In this Letter, we are going to advocate a method to systematically improve the 
self-consistent harmonic approximation (or the Gaussian approximation), which has been 
widely used in condensed matter physics. (See Ref. [6] and references therein.) Indeed 
such a method was applied to the IIB matrix model up to the 3rd order in Ref. [7]. This 
provided the first analytical evidence for the dynamical generation of 4d space-time, which 
is both surprising and encouraging since the model describes superstrings in 10 dimensions. 
Stability of this conclusion against higher order corrections has been confirmed recently by 
explicit 5th order calculations [8]. There the method is also interpreted as an improved 
Taylor expansion. The Gaussian approximation (the leading order of the Gaussian ex- 
pansion) has been applied to random matrix models earlier in Ref. [9]. Its application to 
supersymmetric matrix quantum mechanics (including the Matrix Theory) was advocated 
by Kabat and Lifschytz [10], and Refs. [11] succeeded in revealing interesting blackhole 
thermodynamics even at the leading order. 

Here we apply the method to the bosonic version of the IIB matrix model. More 
generally we consider the model [12] obtained from dimensional reduction of SU(iV) Yang- 
Mills theory in D dimensions, where the D = 10 case corresponds to the bosonic IIB matrix 
model. The model can also be regarded as a Hermitian matrix version of the Eguchi-Kawai 
model [13], which was proposed as an equivalent description of large- N gauge theory. This 
equivalence in the case of the present model has been investigated recently [14]. Hence we 
expect that our results would also have some implications to field theoretical applications 
of the method. 

In the present bosonic model, we perform calculations up to the 7th order, which 
allows us to convincingly demonstrate the convergence of the method. Although the model 
is much simpler than the IIB matrix model, it is still nontrivial, and in particular it has not 
been solved exactly so far. Also the conventional perturbative approach cannot be applied 
as in the IIB matrix model, since the action does not contain a quadratic term. However, 
standard Monte Carlo simulation is applicable, and there are some results [15] that can 
be compared with ours. As another advantage of studying this model, we can perform 
analytic calculations leaving D as a free parameter. We observe that the convergence is 
fast in particular for large D. For D > 10, for instance, a reasonable convergence is already 
achieved at the 3rd order, which supports the validity of the aforementioned calculations 
in the IIB matrix model. 

We would like to point out that actually the Gaussian expansion method has a long his- 
tory [16, 17], although its application to matrix models is new. Originally it was developed 
to obtain the energy spectrum in quantum mechanical systems. The most striking aspect 
of the method is that the results converge, although the expansion is not based on any small 
parameter 1 . It is also interesting that the results contain fully nonperturbative effects, 
although the actual calculations are nothing but the familar ones in standard perturbation 

1 There are rigorous proofs of convergence in certain examples [18]. Very recently Ref. [19] has shown 
that complex solutions to the rf(result) /<i(unphysical parameter) = equation (in fact, the solution with 
the largest imaginary part) works extremely well. 
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theory. In this approach one typically encounters a situation in which the approximants 
depend on free parameters. In fact the power of the method comes from the flexibility 
to adjust these parameters depending on the order of the expansion. Common wisdom 
suggests to determine them by 'the principle of minimum sensitivity' [17]. This is usu- 
ally achieved by solving the 'self-consistency equation'. However, we will encounter some 
problems with this strategy. Here we propose a novel prescription based on histograms, 
with which one can naturally obtain the 'best approximation' at each order without such 
problems. We hope that our histogram prescription is useful in other applications as well. 

2. Dimensionally reduced Yang-Mills integrals 

The model we study in this work is defined by the partition function 

d^e" 5 , (2.1) 



S = -\N[3^2tr[A^,A u ] 2 , (2.2) 

where A^ (// = 1, ■■■,D) are N x N Hermitian matrices. The integration measure dA 
is defined by d^4 = \\ a=1 n^=i "^f > where A® are the coefficients in the expansion 
= Y, a A ti Ta with respect to the SU(iV) generators T a normalized as tr (T a T b ) = \5 ab . 
The model can be regarded as the zero-volume limit of SU(iV) Yang-Mills theory in D 
dimensions. As a result of this limit, we can actually absorb (3 in (2.2) by rescaling 

\ ^ /?~ 1/4 ^V Thus f 3 

is not so much a coupling constant as a scale parameter, which 
we set to (3 = 1 without loss of generality. The partition function was conjectured [12] and 
proved [20] to be finite for N > D/(D — 2). Note in particular that the model is ill-defined 
for D < 2. A systematic 1/D expansion was formulated in [15]. In particular the absence 
of SO(-D) breaking was shown to all orders of the 1/D expansion and this conclusion was 
also confirmed by Monte Carlo simulations [15] for various D = 3, 4, 6, ■ ■ ■ , 20. The Gaus- 
sian expansion method is also capable of addressing such an issue [7], but here we assume 
the absence of the SSB in order to focus on the convergence of the method itself. In Ref. 
[21] Polyakov line and Wilson loop were calculated by the (next-)leading Gaussian approx- 
imation, and Monte Carlo results [3] were reproduced qualitatively. An equally successful 
result was obtained in the supersymmetric case [22]. 

3. The Gaussian expansion method 

The Gaussian action which has both SU(iV) and SO(L>) symmetries can be written as 

*o = ^I>(V > (3-D 

,u=l 

where v is a positive real parameter. (The constant factor is introduced for conve- 
nience.) Then we rewrite the partition function (2.1) as 

Z = Z (e-( s - 5 °))o , (3.2) 
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Z = J dAe~ s ° , (3.3) 

where ( ■ )o is a VEV with respect to the partition function Zq. From (3.2) it follows that 
the free energy F = — In Z can be expanded as 

oo 

F = Y J ?k ; F = -lnZ , (3.4) 

fc=0 

F k = -^r(( S ~ S o) k )c,o (for k > 1) , (3.5) 

where the suffix 'C in ( • )c,o means that the connected part is taken. Each term Fk 
in the above expansion can be calculated conveniently by using Feynman diagrams. In 
actual calculations we have to truncate the infinite series (3.4) at some finite order. In Ref. 
[8], Schwinger-Dyson equations were used to reduce the number of diagrams considerably. 
Using this technique, we were able to proceed up to the 7th order in the present model 
with reasonable efforts (e.g., we evaluated 21 eight-loop diagrams.). 
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Figure 1: The truncated free energy density / is plotted as a function of v for D = 10. Each curve 
corresponds to order 1, 2, ■ ■ •, 7. Formation of a plateau is clearly seen. 

Once we truncate the expansion (3.4), the calculated free energy depends on the pa- 
rameter v introduced in the Gaussian action (3.1). In the large- N limit 2 we compute the 
'free energy density' defined by 



/ 



lim 



D(N 2 - 1) 




The second term is subtracted in order to make the quantity finite. (In particular / = 
for D = oo.) In Fig. 1 we plot the result at each order as a function of v for D 



(3.6) 

-1/4 
= 10. 



2 This simply amounts to considering only the planar diagrams. Although it reduces our effort consider- 
ably, the method itself is applicable to arbitrary finite N. 
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Formation of a plateau is clearly seen. We may understand this phenomenon as a reflection 
of the fact that the formal expansion (3.4), if not truncated, should not depend on v. 
According to this interpretation, the height of the plateau observed for the truncated free 
energy is expected to provide a good approximation. This is the philosophy behind 'the 
principle of minimum sensitivity', which is known in a more general context [17]. 
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Figure 2: The zoom up of Fig. 1 for order 7. The arrows indicate the position of the extrema, 
each of which corresponds to a solution of the 'self-consistency equation'. 



Next we need to specify a prescription to obtain a 
concrete value which approximates the free energy at 
each order. In Refs. [7, 8] the free parameters in the 
Gaussian action were determined such that the trun- 
cated free energy becomes stationary. This amounts 
to solving the 'self-consistency equation' 
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Table 1: The number of solutions to 
the 'self-consistency equation' (3.7) for 
each D. 



In Table 1 we list the number of solutions to the 'self- 
consistency equation' for each D, which already re- 
veals some problems. For example, for D = 10 we 
find that the number of solutions increases as we go 
to higher orders, and we have to decide which one to 
choose. A similar ambiguity was encountered in the IIB matrix model [7, 8], where the 
solution that gives the smallest free energy has been chosen. While this prescription seems 
reasonable at relatively small orders, we have problems at higher orders, in particular when 
we are discussing the convergence. In Fig. 2 we zoom up the plateau seen at the order 7 
in Fig. 1. We observe some oscillations, whose amplitude becomes larger towards the left 
edge of the plateau. If we used the prescription described above in the present situation, 
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we would clearly underestimate the height of the plateau. This will be even more prob- 
lematic at higher orders since the dip at the left edge seems to become deeper as the order 
increases. It may also happen in general that the truncated free energy acquires a global 
minimum which has nothing to do with the plateau formation. 

From Table 1 one also finds that there are cases where the 'self-consistency equation' 
has no solutions. (This was encountered in the study of the IIB matrix model at the orders 
2 and 4 [7, 8].) For instance there is no solution at the orders 2,4,6 for D = 4. As we show 
in the Appendix, a well-defined plateau is developing even in this case. The question is 
how to extract the height of the plateau in such cases. 

4. The histogram prescription 

All these problems arise from trying to determine the free parameters by local information 
such as 'self-consistency equations'. In fact it is more natural to extract the height of 
the plateau directly. As a first step, let us make a histogram of the truncated free energy 
calculated at points which are uniformly distributed in the parameter space of the Gaussian 
action. Fig. 3 shows the histogram for D = 10 at order 7. Whenever we make histograms 
in what follows, the range of v is restricted to 0.01 < v < 1.4, and the vertical axis is 
normalized such that the integration gives unity. In Fig. 3 the bin size is chosen to be 
very small, and as a result we obtain many peaks, each of which corresponds to a solution 
of the 'self-consistency equation' (Actually the highest peak contains the two solutions 
which lie at the right edge of the plateau in Fig. 2.). Hence we are back in the previous 
situation, and we cannot reasonably extract the height of the plateau. In order to identify 
the plateau, we clearly need some kind of coarse graining. Let us note here that the bin size 
of the histogram represents the precision at which we search for the plateau. Therefore, 
we increase the bin size gradually until a single peak becomes 'dominant'. As a criterion 
of the dominance, we require the highest bin to be more than twice as high as the second 
highest one. We further require the dominance to persist for larger bin size. The resulting 
histogram is shown in Fig. 4. 

Fig. 5 shows the histograms obtained in this way at each order for D = 10. The 
peak becomes sharper and sharper as we go to higher order. Thus the histogram serves 
as a clear indicator of the plateau formation. Moreover, now we can obtain the 'best 
approximation' at each order from the position of the peak, while the bin size (chosen 
by the above criterion) roughly represents the theoretical uncertainty of the method. In 
the Appendix, we present a similar analysis for D = 4. In this case the 'self-consistency 
equation' has no solutions at orders 2,4,6 as mentioned above. However, the histogram 
prescription is as successful as in D = 10. 

In Fig. 6, we plot the approximated value of the free energy for various D obtained at 
each order by our histogram prescription. The 'error bars' represent a rough estimate for 
the theoretical uncertainty explained above 3 . (Remark: Unlike usual error bars, they do 

3 We did not put error bars to the order 1 results. This is because our criterion for the dominance of 
a single bin is satisfied even for an infinitesimal bin size. Thus at the order f, the histogram prescription 
does nothing more than just solving the 'self-consistency equation'. 
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Figure 3: The histogram of the truncated free energy density / for D = 10 at order 7. The bin 
size is chosen to be very small (0.00002). 
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Figure 4: The histogram of the truncated free energy density / for D — 10 at order 7. The bin 
size (0.00035) is chosen such that the highest bin is more than twice as high as the second highest 
one. 



not represent an estimate for possible discrepancies from the correct value. For the latter, 
one should also consider the changes of the results as the order is increased.) We observe 
a clear convergence for each D. The convergence becomes fast in particular for large D. 
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Figure 5: The histogram of the truncated free energy density / for D = 10 at order 1, 2, • • •, 7. 
The bin size is chosen carefully as described in the text at each order. The growth of the peak is 
clearly observed. 
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Figure 6: Free energy density / obtained by the histogram prescription at order 1, 2, ■ ■ ■, 7 for 
various D. The horizontal line represents the exact result / = —1/4 for D = oo. 

5. Observables 

Similarly to the free energy (3.4), the expectation value of an operator O is expanded as 

oo 

(0) = J2°k ; O ^(0)o, (5.1) 

(5.2) 



A;=0 
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As a fundamental object in this model [15], we consider 




(5.3) 



The normalization is chosen such that (O) = 1 for D = oo in the large- N limit [15]. We 
calculate (O) by the series expansion (5.1) truncated at some order. In Fig. 7 we plot the 
result at each order as a function of v for D = 10. 

In Refs. [7, 8], the truncated expectation value was evaluated at the same value of 
v as the one used to evaluate the free energy. Note, however, that the region of the 
parameter v where the plateau develops for the expectation value is slightly shifted from 
the corresponding region for the free energy. It is therefore safer to extract an approximated 
value for the observable by repeating the same procedure as we did for the free energy. 

Fig. 8 shows the results of the histogram prescription at each order. The growth of the 
peak is clearly observed, which confirms the plateau formation. Fig. 9 summarizes the final 
results for the observables, where the 'error bars' again represent the theoretical uncertainty. 
Comparison with the Monte Carlo data [15] confirms that our results do converge to the 
correct values. For D > 10 in particular, the convergence is already achieved at order 3. 

It was previously noted that the leading Gaussian approximation becomes 'exact' at 
D = oo [21]. The precise statement is that if we evaluate the free energy density and 
the observable (truncated at the 1st order) at the value of v which solves the 1st order 
'self-consistency equation' for the free energy, they agree with the exact results at D = oo. 
Although this can be naturally understood from the viewpoint of the 1/D expansion [15], it 
looks rather accidental in the present framework. If we evaluate the truncated observable at 
v which solves the 'self-consistency equation' for the observable, we obtain 1.08866 instead 
of 1. Also there is no plateau yet at order 1. Figs. 6 and 9 show that the convergence of 
the method as it stands is achieved at order 3 (but not yet at order 1) even for D = oo. 

In this regard, it should be emphasized that the Gaussian expansion is different from 
the 1/D expansion. Each term in the expansions for the free energy density and the 
observable is of order 1 with respect to 1/D for generic v. Still the height of the plateau 
converges due to nontrivial cancellations among diagrams. Note also that the Gaussian 
expansion can be applied to the supersymmetric case (including the IIB matrix model), 
while the 1/D expansion cannot be formulated there. 

6. Summary and Discussions 

The most important conclusion of our study is that the Gaussian expansion method con- 
verges in the dimensionally reduced Yang-Mills integrals. In particular the convergence 
is achieved already at order 3 for D = 10, which corresponds to the bosonic IIB matrix 
model. This suggests the particular usefulness of the method in studying the IIB matrix 
model. In demonstrating the convergence we encountered problems associated with the 
commonly used prescription based on 'self-consistency equations'. We solved these prob- 
lems by extracting the height of the plateau directly using a novel prescription based on 
histograms. 
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Figure 7: The truncated expectation value (O) is plotted as a function of v for D = 10. Each 
curve corresponds to order 1, 2, • • ■, 7. Formation of a plateau is clearly seen. 
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Figure 8: The histogram of the truncated expectation value (O) for D = 10 at order 1, 2, ■ ■ ■, 7. 
The bin size is chosen carefully as described in the text at each order. The growth of the peak is 
clearly observed. 



The histogram prescription is expected to be particularly useful in the case where the 
parameter space of the Gaussian action becomes multi-dimensional. If there are more than 
two parameters, it is already difficult to visualize the plateau formation as we did in Figs. 
1, 7 and 10. When we study the SSB of SO(10) symmetry in the IIB matrix model, we 
have to introduce 10 real and 120 complex parameters [7]. Here we will need some kind of 
important sampling technique like Monte Carlo simulations in order to explore such a huge 
parameter space. In the previous works [7, 8] based on 'self-consistency equations', there 
were problems like the ones we saw in this work. At order 2 and 4 there were no solution to 
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Figure 9: The expectation value (O) obtained by the histogram prescription at order 1, 2, ■ ■ -, 
7 for various D. The horizontal lines represent the Monte Carlo results for each D < oo and the 
exact result (O) — 1 for D = oo. 



the 'self-consistency equations', and in the other cases the solution was not unique. We are 
currently reinvestigating the IIB matrix model by using the histogram prescription with 
the above extension [23]. 

Finally we would like to emphasize that the dimensionally reduced model we studied 
in this paper is a system of infinitely many degrees of freedom, since we have taken the 
size of the matrices to be infinite. In this regard, let us recall that the model has also 
connections to SU(iV) gauge theories in D-dimensional space time in the large- N limit [13]. 
The dimensional reduction amounts to reducing the base space of the SU(iV) gauge theories 
to a point, but the original space-time degrees of freedom are still somehow encoded in the 
internal degrees of freedom. This is possible precisely because there are infinitely many 
internal degrees of freedom after taking the large- N limit. From this point of view the 
convergence of the Gaussian expansion method in the present model is interesting since 
the examples for which the convergence is shown in the literature have been restricted to 
simple quantum-mechanical systems. Although higher order calculations in field-theoretical 
applications would be technically more involved due to the necessity of the renormalization 
procedure (See, however, Ref. [24]), we hope that our results open up a new perspective in 
that direction as well. 
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A. The situation in D = 4 : plateau without extrema 

In this Appendix, we present our results for D = 4. As we see from Fig. 10, a well-defined 
plateau is clearly developing. However, if we zoom up the curves for order 2,4,6 (Fig. 
11), we find that there are no extrema, in agreement with the absence of solutions to the 
'self-consistency equation' (See Table 1). We will see that the histogram prescription works 
even in this case. 

Figs. 12 and 13 show the histograms at order 6 before and after 'coarse-graining'. As 
we increase the bin size gradually, a single bin becomes dominant (See Section 4 for the 
precise criterion.). Thus the histogram prescription allows us to obtain explicit results at 
each order including orders 2,4,6, where the 'self-consistency equation' has no solutions. 
Fig. 14 shows the results at orders 1, ■ ■ ■ , 7. The growth of the peak is clearly observed, 
which confirms the plateau formation. Moreover, the results obtained in this way show a 
clear convergence as we have seen in Fig. 6. 
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Figure 10: The truncated free energy density / is plotted as a function of v for D = 4. Each curve 
corresponds to order 1, 2, • • ■, 7. Formation of a plateau is clearly seen. 
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Figure 13: The histogram of the truncated free energy density / for D = 4 at order 6. The bin 
size (0.002) is chosen such that the highest bin is more than twice as high as the second highest 
one. 
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Figure 14: The histogram of the truncated free energy density / for D = 4 at order 1, 2, ■ ■ •, 7. 
The bin size is chosen carefully as described in the text at each order. The growth of the peak is 
clearly observed. 
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